function [ber_time_div, diversity_order, fig_handle] = time_diversity(EbN0_dB, num_bits, time_slots)
% 时间分集性能分析
% 输入:
%   EbN0_dB - 比特信噪比范围 (dB)
%   num_bits - 仿真比特数
%   time_slots - 时间分集时隙数
% 输出:
%   ber_time_div - 时间分集误比特率矩阵
%   diversity_order - 分集阶数
%   fig_handle - 图形句柄

    EbN0_linear = 10.^(EbN0_dB/10);
    ber_time_div = zeros(length(EbN0_dB), length(time_slots));
    
    % 添加路径
    addpath('../Common');
    colors = color_definitions();
    
    fprintf('=== 时间分集性能分析 ===\n');
    
    for div_idx = 1:length(time_slots)
        L = time_slots(div_idx); % 分集阶数
        fprintf('仿真时间分集阶数 L = %d...\n', L);
        
        for eb_idx = 1:length(EbN0_dB)
            fprintf('  Eb/N0 = %d dB...\n', EbN0_dB(eb_idx));
            
            % 生成随机比特
            tx_bits = randi([0, 1], num_bits, 1);
            
            % BPSK调制
            tx_symbols = 2 * tx_bits - 1;
            
            % 时间分集：重复发送L次
            rx_symbols_combined = zeros(num_bits, 1);
            
            for l = 1:L
                % 瑞利衰落信道 (每个时隙独立)
                rayleigh_fading = sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
                
                % 添加AWGN噪声
                noise_power = 1 / (2 * L * EbN0_linear(eb_idx)); % 总功率分配
                noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
                
                % 接收信号
                rx_symbols = tx_symbols .* rayleigh_fading + noise;
                
                % MRC合并 (最大比合并)
                weights = conj(rayleigh_fading) ./ (abs(rayleigh_fading).^2 + noise_power);
                rx_symbols_combined = rx_symbols_combined + weights .* rx_symbols;
            end
            
            % 判决
            rx_bits = (real(rx_symbols_combined) > 0);
            
            % 计算误比特率
            ber_time_div(eb_idx, div_idx) = mean(tx_bits ~= rx_bits);
        end
    end
    
    % 计算分集阶数
    diversity_order = zeros(1, length(time_slots));
    for div_idx = 1:length(time_slots)
        % 使用高SNR区域的斜率估计分集阶数
        high_snr_idx = EbN0_dB > 10; % SNR > 10 dB
        if sum(high_snr_idx) > 2
            slope = diff(log10(ber_time_div(high_snr_idx, div_idx))) ./ diff(EbN0_dB(high_snr_idx));
            diversity_order(div_idx) = -mean(slope(:)); % 负斜率，确保标量
        else
            diversity_order(div_idx) = time_slots(div_idx); % 理论值
        end
    end
    
    % 绘制结果
    fig_handle = figure('Name', '时间分集性能', 'Position', [200, 200, 1200, 800]);
    
    % BER性能
    subplot(2,3,1);
    for div_idx = 1:length(time_slots)
        semilogy(EbN0_dB, ber_time_div(:, div_idx), ...
                 'LineWidth', 2, 'Color', colors(div_idx,:));
        hold on;
        
        % 理论曲线 (L阶分集)
        theory_ber = nchoosek(2*time_slots(div_idx)-1, time_slots(div_idx)) * ...
                     (0.5 * (1 - sqrt(EbN0_linear ./ (1 + EbN0_linear)))).^time_slots(div_idx);
        semilogy(EbN0_dB, theory_ber, 'LineWidth', 1, 'Color', colors(div_idx,:), 'LineStyle', '--');
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('时间分集BER性能');
    legend([arrayfun(@(x) sprintf('L=%d', x), time_slots, 'UniformOutput', false), ...
           arrayfun(@(x) sprintf('L=%d(理论)', x), time_slots, 'UniformOutput', false)]);
    
    % 分集阶数验证
    subplot(2,3,2);
    plot(time_slots, diversity_order, 'o-', 'LineWidth', 2, 'Color', colors(1,:));
    hold on;
    plot(time_slots, time_slots, 'k--', 'LineWidth', 2); % 理论直线
    grid on;
    xlabel('时间分集阶数 L');
    ylabel('实测分集阶数');
    title('分集阶数验证');
    legend('实测值', '理论值');
    
    % 分集增益
    subplot(2,3,3);
    ber_ref = ber_time_div(:, 1); % 无分集作为参考
    for div_idx = 2:length(time_slots)
        div_gain = ber_ref ./ ber_time_div(:, div_idx);
        semilogy(EbN0_dB, div_gain, 'LineWidth', 2, 'Color', colors(div_idx,:));
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('分集增益');
    title('分集增益 vs SNR');
    legend(arrayfun(@(x) sprintf('L=%d vs L=1', x), time_slots(2:end), 'UniformOutput', false));
    
    % 固定BER下的SNR增益
    subplot(2,3,4);
    target_ber = 1e-3;
    snr_gains = zeros(1, length(time_slots));
    
    for div_idx = 1:length(time_slots)
        % 插值找到目标BER对应的SNR
        ber_vec = ber_time_div(:, div_idx);
        % 确保BER值是唯一的（避免重复值）
        [unique_ber, unique_idx] = unique(ber_vec);
        unique_ebno = EbN0_dB(unique_idx);
        
        if min(ber_vec) <= target_ber && max(ber_vec) >= target_ber
            snr_gains(div_idx) = interp1(unique_ber, unique_ebno, target_ber);
        else
            snr_gains(div_idx) = NaN;
        end
    end
    
    % 计算相对于无分集的增益
    if ~isnan(snr_gains(1))
        snr_gains = snr_gains - snr_gains(1);
        plot(time_slots, snr_gains, 'o-', 'LineWidth', 2, 'Color', colors(1,:));
        hold on;
        grid on;
        xlabel('时间分集阶数 L');
        ylabel('SNR增益 (dB)');
        title(sprintf('SNR增益 @ BER=%.0e', target_ber));
    end
    
    % 信道时间相关性影响
    subplot(2,3,5);
    correlation_values = [0, 0.5, 0.9, 0.99];
    L = 4; % 固定分集阶数
    ber_correlated = zeros(length(EbN0_dB), length(correlation_values));
    
    for corr_idx = 1:length(correlation_values)
        rho = correlation_values(corr_idx);
        fprintf('仿真时间相关性 rho = %.2f...\n', rho);
        
        for eb_idx = 1:length(EbN0_dB)
            tx_bits = randi([0, 1], num_bits, 1);
            tx_symbols = 2 * tx_bits - 1;
            
            % 生成相关信道
            channel_states = zeros(num_bits, L);
            channel_states(:, 1) = sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            
            for l = 2:L
                innovation = sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
                channel_states(:, l) = rho * channel_states(:, l-1) + sqrt(1-rho^2) * innovation;
            end
            
            % 合并
            rx_combined = zeros(num_bits, 1);
            for l = 1:L
                noise_power = 1 / (2 * L * EbN0_linear(eb_idx));
                noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
                rx_symbols = tx_symbols .* channel_states(:, l) + noise;
                
                weights = conj(channel_states(:, l)) ./ (abs(channel_states(:, l)).^2 + noise_power);
                rx_combined = rx_combined + weights .* rx_symbols;
            end
            
            rx_bits = (real(rx_combined) > 0);
            ber_correlated(eb_idx, corr_idx) = mean(tx_bits ~= rx_bits);
        end
    end
    
    for corr_idx = 1:length(correlation_values)
        semilogy(EbN0_dB, ber_correlated(:, corr_idx), ...
                 'LineWidth', 2, 'Color', colors(corr_idx,:));
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('时间相关性影响 (L=4)');
    legend(arrayfun(@(x) sprintf('ρ=%.2f', x), correlation_values, 'UniformOutput', false));
    
    % 功率分配策略
    subplot(2,3,6);
    power_allocations = {'等功率', '注水算法', '自适应'};
    ber_power = zeros(length(EbN0_dB), length(power_allocations));
    L = 4;
    
    for alloc_idx = 1:length(power_allocations)
        fprintf('仿真功率分配: %s...\n', power_allocations{alloc_idx});
        
        for eb_idx = 1:length(EbN0_dB)
            tx_bits = randi([0, 1], num_bits, 1);
            tx_symbols = 2 * tx_bits - 1;
            
            % 生成信道
            channels = zeros(num_bits, L);
            for l = 1:L
                channels(:, l) = sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            end
            
            % 功率分配
            if alloc_idx == 1 % 等功率
                powers = ones(1, L) / L;
            elseif alloc_idx == 2 % 注水算法
                gains = mean(abs(channels).^2, 1); % 平均信道增益
                [powers, ~] = waterfilling(gains, 1, 1/(2*EbN0_linear(eb_idx)));
            else % 自适应 (基于瞬时信道)
                powers = zeros(1, L);
                for bit_idx = 1:num_bits
                    gains = abs(channels(bit_idx, :)).^2;
                    [powers_inst, ~] = waterfilling(gains, 1, 1/(2*EbN0_linear(eb_idx)));
                    powers = powers + powers_inst;
                end
                powers = powers / num_bits;
            end
            
            % 发送和合并
            rx_combined = zeros(num_bits, 1);
            for l = 1:L
                noise_power = 1 / (2 * powers(l) * EbN0_linear(eb_idx));
                noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
                rx_symbols = tx_symbols .* channels(:, l) + noise;
                
                weights = conj(channels(:, l)) .* sqrt(powers(l)) ./ (abs(channels(:, l)).^2 + noise_power);
                rx_combined = rx_combined + weights .* rx_symbols;
            end
            
            rx_bits = (real(rx_combined) > 0);
            ber_power(eb_idx, alloc_idx) = mean(tx_bits ~= rx_bits);
        end
    end
    
    for alloc_idx = 1:length(power_allocations)
        semilogy(EbN0_dB, ber_power(:, alloc_idx), ...
                 'LineWidth', 2, 'Color', colors(alloc_idx,:));
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('功率分配策略 (L=4)');
    legend(power_allocations);
    
    % 保存图片
    saveas(fig_handle, 'results/time_diversity_analysis.png');
end

function [powers, capacity] = waterfilling(gains, total_power, noise_power)
% 注水算法功率分配
    mu = 1; % 初始水位
    powers = zeros(1, length(gains));
    
    % 迭代找到最优水位
    for iter = 1:100
        powers = max(0, mu - noise_power ./ gains);
        if abs(sum(powers) - total_power) < 1e-6
            break;
        end
        mu = mu + (total_power - sum(powers)) / sum(powers > 0);
    end
    
    % 计算容量
    capacity = sum(log2(1 + powers .* gains / noise_power));
end